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Abstract We present filling as a new type of spatial subdivision problem that 
is related to covering and packing. Filling addresses the optimal placement of 
overlapping objects lying entirely inside an arbitrary shape so as to cover the 
most interior volmxie. In n-dimensional space, if the objects are polydisperse 
n-balls, we show that solutions correspond to sets of maximal n-balls and the 
solution space can reduced to the medial axis of a shape. We examine the 
structure of the solution space in two dimensions. For the filling of polygons, 
we provide detailed descriptions of a heuristic and a genetic algorithm for 
finding solutions of maximal discs. We also consider the properties of ideal 
distributions of N discs in polygons as — ^ oo. 

Keywords Polygons • Filling • Packing • Covering 

PACS 47.57.J- • 02.70.-c • 47.57.Bc • 88.80.ht 

Mathematics Subject Classification (2010) 05C70 • 52C15 • 52C17 
1 Introduction 

First introduced in reference we define filling as the problem of packing 
overlapping objects inside of a defined shape such as to optimally cover the 
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interior volume without extending beyond the boundary of the shape. We are 
primarily interested in the optimal filling of an n-dimensional shape with a 
well-defined n — 1 surface with n-dimensional polydisperse balls. 

The filling problem can be expressed by the following two questions: 

Problem 1 Given a compact region G (having non-empty interior and no 
holes in the interior) and a fixed positive integer TV, how can N balls of varying 
radii be placed completely interior to G so as to maximize the total volume 
covered? Overlaps of the N balls are permitted. 

Problem 2 In general, for each fixed shape G, what is the best strategy for 
maximizing the fraction of volume covered by the minimum number of balls 
of varying radii in G? 

In the deceptively simple problem of determining the optimal set of balls 
to fill an arbitrary shape we find a surprisingly rich problem with many open 
questions. In this paper we address the above two questions. In Section [2] 
we define the filling problem and the basic terminology necessary to discuss 
the problem. We describe the mathematical structure of the filling solution 
space and show that it can be reduced from dimension n -I- 1 to dimension 
n—1. We characterize the forms of degeneracy possible in filling solutions and 
demonstrate how each individual ball contributes to the filling of a shape. In 
Section[3j we first restrict the filling problem to planar shapes. We then identify 
features of the solution space that affect the properties of an optimal solution 
and the methods for finding solutions. We introduce the concept of neighbors 
and show how a fixed point in the solution space divides the solution space 
into independent spaces. We show how special points in the solution space are 
found in many solutions. We then restrict the problem further to polygons, 
where the structure of the solution space can be reduced to a small number of 
cases. We numerically explore the solution space of a simple construction of 
three discs, and show the solution space is complex with many local maxima 
and topologically diverse configurations. In Section|4]we detail two algorithms 
for generating filling solutions for polygons, a genetic algorithm and a heuristic 
algorithm. The genetic algorithm utilizes a minimal set of assumptions about 
the solution space. The heuristic algorithm exploits the known structure of the 
solution space and also relies on several conjectures. We discuss the relative 
efficiency of the heuristic algorithm in searching the solution space and the 
good correspondence between the two algorithms. In Section [5] we find the 
distribution of discs in a polygon at the continuum limit, or as iV — >■ oo. The 
derived analytical expressions may be used to approximate solutions for finite 
but large N. We derive an expression for the fractional allocation of discs over 
the medial axis branches an arbitrary polygon and an exact expression for a 
triangle as iV — )■ oo. In Section [6j we provide concluding remarks. In Section 
[Sj we include a glossary of terms defined in this paper. 
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2 General Properties of the Medial Axis of G and Filling Solutions 

2.1 Definitions and Theorems 

Let G be a compact (closed and bounded), simply-connected n-dimensional 
region with a non-empty interior. Let S be the boundary oi G: S = 5 G. 
As in reference [5], we restrict S to have a tangent and curvature defined 
everywhere but at a finite number of points. At these points, sided curvature, 
i.e. a limit performed only on one side of the point, exists from any direction 
along the boundary. For simplicity, we do not consider G with holes, nor G 
with a boundary that abuts itself. 

Definition 1 Let i? at be a set containing N balls Di that are completely 
contained in G. Each Di has a radius and center Xj. i?jv is a filling solution 
of G. Let (f>{RN, G) be the fraction of G that is covered by Rn- 4> is the measure 
of the filling G by the set Rm- 

The measure (j) is equal to the volume of the union of balls of Rn divided 
by the volume of G, thus (/> < 1 by definition, (j) is equal to unity for N < oo 
only if G is equivalent to a finite number of overlapping balls. The space of all 
Rn is of dimension n + 1. 

Definition 2 If, VA e Rn, HRn ~ {A},G) < (j){RN,G), then the set is 
all-filling. In other words, each ball in Rn uniquely fills a non-zero volume of 
G. 

Definition 3 A set Rn is an optimal filling solution of G if there is no other 
set R'j^ that satisfies (j){R'^,G) > (I){Rn,G). 

The function (p can be defined over the space of all _Rjv for a shape G and 
fixed positive N. Our objective is to find the set Rn with the maximum value 
of (j)- We will now prove that the solution space can be restricted to sets of Rn 
containing only maximal n-balls, a space defined by the medial axis of G and 
its associated radius function. 

The medial axis of an object, originally defined by Blum in reference [3], 
and also known as the topological or medial skeleton, is the set of all points 
having more than one closest point on the object's boundary. We use the 
notation M(G) for the medial axis of G. The medial axis is a reduction of an 
n-dimensional shape into an n — 1-dimensional space, the locus of centers of 
the maximal n-balls. A maximal n-ball is defined as a ball that is tangent to 
the boundary at two or more points. It is also a ball contained completely in 
G that is not a proper subset of any other ball also contained in G. A shape 
is the logical union of all its maximal n-balls. The radius function associated 
with M(G) is a continuous, non-negative function defined at each point of 
M(G) as the radius of the maximal n-ball centered at that point. The medial 
axis and the radius function together are a complete shape descriptor [2] and 
can be used to reconstruct the shape. 
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Theorem 1 ForG, there exists optimal filling solutions Rn that contain only 
maximal balls. 

Proof From any filling solution set that has a ball that is not on the medial 
axis, we can construct a solution set that contains only balls on the medial axis. 
Assume we have a solution set Rn that contains a ball D that is not tangent 
to S at any point. That ball is completely contained inside a concentric ball 
that is tangent to at least one point of S. And that ball is completely contained 
inside a larger cotangent ball that is also tangent to a second point of S. This 
last ball D' has its center on some part of the medial axis by construction and 
is thus a maximal ball. Let R'^ be the set of balls where D is replaced by D' . 
It must be that 4){R'^, G) > (I>{Rn, G). So if Rn is an optimal filling solution, 
then so is R'n- 

While Theorem [T] implies that filling solutions can be restricted to sets of 
maximal balls, it does not follow that optimal solutions must be composed of 
maximal balls for all shapes. Shapes that have boundaries with concave points 
of infinite curvature can have optimal filling solutions with non-maximal balls. 
Figure [ija) shows such a shape. The outer boundary of the shape in Figure 
[ija) is defined by four circular arcs. The dashed (green online) line is the 
medial axis of this shape. If the four discs at the extreme points of the medial 
axis have been placed, then there is no need for the final disc placed inside the 
shape to be a maximal disc. 

Even when restricted to optimal solutions of maximal balls, the entire 
medial axis need not be occupied as N ^ oo. Figure [ijb) is an example of a 
concave shape with two concave points. The portion of the medial axis between 
the two circle centers (red online) need not be occupied by disc centers to fill 
the shape as — >■ cx). 

Theorem 2 If S contains no concave points of infinite curvature, then op- 
timal fillings Rn composed of maximal balls are also all-filling. Only filling 
solutions of maximal balls can be optimal. The entire medial axis is occupied 
for optimal filling solutions as N oo. 

Proof Assume there is an all-filling filling solution Rn for a shape G that 
contains a ball D that is not a maximal ball. The operations from Theorem [l] 
are used to construct a maximal ball D' on the medial axis from the ball D. 
The disc £>', by construction, is tangent to the boundary of S in at least one 
location that D was not. Assuming that S has no concave points of infinite 
curvature (e.g a reflex vertex for a polytope), then the point of tangency is 
smooth and there is a small region around the point of tangency that D' covers 
that D did not (Figure [2]) . That region can only already have been covered by 
a ball in Rn if a ball contained in Rn of equal or larger radius, was tangent 
to S the same point. But since disc D' , and therefore D, would be completely 
contained in that ball, then Rn would not be all-filling. So (j){R'N, S) is strictly 
greater than <f>{RN, S). If each point on the boundary S is tangent to a maximal 
ball at a smooth point, then there is one maximal ball tangent to S at that 
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point, so each point of S maps to a one and only one point on the medial axis. 
Also, a unique infinitesimal volume is covered by that maximal ball, so the 
entire medial axis must be occupied for the optimal fillings as — ^ oo. 




(a) (b) 

Fig. 1 (a)Tho construction has optimal solutions without maximal balls. The center of 
the red disc need not be on the medial axis (dashed green) for the shape to be completely 
covered, (b) The construction need not have all of its medial axis (dashed green) filled. A 
disc added to the red portion fills no additional area. 



Fig. 2 If the point of tangency is smooth in any direction, the largest ball tangent to the 
point covers more volume than any smaller ball tangent to the point. 



Theorem [T] shows that to construct filling solutions, the search space can 
be restricted to the space of maximal balls. Theorem [2] shows that for G with 
S without concave points of infinite curvature, optimal filling solutions consist 
only of maximal balls. Searching for optimal fillings has been reduced from 
finding points in an n + 1 dimensional space (disc center position and radius) 
to finding points on an n — 1 dimensional surface, or a problem of dimension 
n - 1. 

In practice, this surface is better described as a set of bounded connected 
surfaces. For example, a planar shape has a medial axis that forms a planar 
graph, a connected set of curves that meet at points. In three dimensions, a 
medial axis is composed of sheets, seams, and j unctions [HinilS]. 

In the following sections we shall implicitly assume all _Rjv solutions being 
discussed are all-filling. 

For an arbitrary shape G, there is no guarantee that optimal filling solu- 
tions must be unique. For example. Figure [3] shows how a long rectangle can 
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(b) 

Fig. 3 Two examples of degenerate solutions. For (a) a single maximal ball can be placed in 
infinitely many locations in a rectangle. For (b), the symmetrical triangle, the asymmetrical 
solution can be reflected to create a degenerate solution. 



have an infinite number of A'' = 1 solutions. All discs added to the middle of 
the rectangle have the same diameter. When sufficiently many discs have been 
added to the rectangle so that the discs must overlap, this form of degeneracy 
disappears. Degenerate solutions also occur in symmetrical shapes with asym- 
metrical optimal filling solutions, such as the N = 2 solution shown for the 
triangle in Figure |3]. This form of degeneracy does not disappear as iV oo. 
It is also possible for a shape with no symmetries to have two distinct optimal 
filling solutions with the same filling. The likelihood of two distinctly different 
sets of discs both to be optimal and have the same is unlikely, and thus this 
form of degeneracy is also likely to be extremely rare. 

2.2 Filling contribution of a single ball 

The contribution of a single ball to the measure (j) is equal to the volume of 
the ball offset by its fractional share of the volume of any overlap with other 
balls. More explicitly, the contribution of a single ball is the volume it uniquely 
covers plus 1/i of the volume it shares with exactly i other balls. Let V^{Dk) 
be the domain that a ball Dk shares with exactly i other balls, including itself. 
Let Mv{V) be the measure of the volume of a domain. Note that V-[{Dk) is 
the domain uniquely covered by the ball Dk- The contribution C{Dk) is, 
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C{D,)=Y,-^Mv(V:{D,)) (1) 

and, 

1 ^ 

In general, if the overlap between balls P and Q is completely contained 
inside of the overlap between balls Q and R, then locally adding, removing, 
or displacing ball P cannot uncover any of that overlap volume. So if ball P 
is added, removed, or locally displaced, the contribution of other balls may 
change but the only term that changes in the summed contributions of all the 
balls is {V(iP)). 



3 Planar Shapes and Polygons 

We now restrict the problem to planar shapes, G, whose medial axes are the 
locus of the centers of maximal discs. Various algorithms exist to compute the 
medial axis of simple polygons (polygons that are not self-intersecting) and 
planar regions bounded by line segments, circular arcs, and general nonuniform 
rational B-spline [ZIEIE]. 

We define the terminology and review the properties for an M{G) of a 
planar shape introduced by Blum and Nagel [5]. M{G) consists of connected 
subsets of points that form a 1-D planar graph. Most points of M{G) are 
normal points, whose maximal disc is in contact with the boundary at two 
separate but contiguous sets of points. M{G) also contains a finite number of 
branch points, each of which has a maximal disc in contact with the boundary 
at three or more separate but contiguous sets of points, and a finite number 
of end points, whose maximal disc is in contact with the boundary at only one 
contiguous set of points. For all but a finite number of discs, the contiguous set 
of points is a single point of contact. The contact point consists of more than 
just a single point if the radius of curvature the disc and boundary are the 
same. As long as G has no holes, then the graph M{G) forms a tree with no 
loops. M{G) can be divided into sets of contiguous normal points bounded by 
branch or end points, such that the division is unique, disjoint, and complete. 
Sets of contiguous normal points shall be referred to as a branch. The boundary 
of S can be divided into parts associated with each branch by the intersection 
of S with the set of maximal discs defined over a branch. This division of S is 
also unique, disjoint, and complete [1]. The shape and radius function of any 
branch of M{G) is determined by the parents of the branch, that is, the two 
contiguous sections of S associated with the given branch. Given a sequence 
of maximal discs on the branch, their respective points of contact with S and 
their centers on the branch are traversed in the same order (assuming the 
branch is traversed in the correct direction relative to S). 
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Definition 4 Given a filling solution Rn, and D G Rn for a AI{G) planar 
graph, we define the neighbors of a maximal disc D, as the maximal discs 
whose centers are the closest along paths in M (G) originating at the center of 
D. 

In other words, if there is a path in M{G) that connects the center of disc D 
to the center of disc D' without traversing another disc center, then disc D 
and disc D' are neighbors (Figure|4|. For G with no holes, where M{G) has no 
loops, there is only one path connecting any two disc centers. When a branch 
of the medial axis is populated with many maximal discs, most of the discs 
have exactly two neighbors, the disc to the left and right of them in sequence. 
If a branch point has connectivity n, then, in a densely populated M{G) the 
disc closest to that branch point has n neighbors. A disc can theoretically have 
as many neighbors as M{G) has end points with connectivity 1. 

Theorem 3 Any overlap of disc D with any disc that is not a neighbor of D 
must be contained inside the overlap of D with one of its neighbors. 

Thus to measure A4v {^[{D)) for a disc D, only the position of the neigh- 
bors of D need be accounted for. To show the latter is true, we draw upon the 
properties of a planar medial axis. 




Fig. 4 Each point represents the center of a maximal disc on M{G). The neighbors of the 
disc A are circled. 



Proof Assume that the centers of three maximal discs A, _B, and C, are on 
a path of M{G). Let their centers be denoted by a, 6, and c. The edges of 
discs A, i?, and G are simply the circles with centers at a, 6, and c of the 
same radii as the discs. Assuming M{G) has no loops, then the path is the 
only path in M{G) connecting the center of A to G. Suppose that there is 
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an intersection between maximal disc A and C as constructed in Fig. [Sj Two 
circles intersect in a region shaped like an asymmetric lens. Label the two 
points where the edge of disc A and C intersect I far and Idose- Divide the 
plane into four quadrants defined by the line connecting a and c and the line 
connecting I far and Idose- Now construct disc B. Without loss of generality, 
we assume the center h is in the bottom right (Quadrant 4) of the figure. As 
disc B is constructed, there are restrictions on both where its center can be 
placed within Quadrant 4 and on its radius. 

First, the point I far must be contained in the disc B. We observe that 
as one is traversing the boundary S of the shape, the edge of A, B, and C 
must be encountered in the following order: a continuous set of disc A edge 
points, a continuous set of disc B edge points, a continuous set of disc C edge 
points, another continuous set of disc C edge points, another continuous set 
of disc B edge points, and another continuous set of disc A edge points. Other 
continuous sets of other disc edge points may interleave the sets specified, 
however, the specified order of encountering continuous sets of points of A, B, 
and C must still be followed. If the radius function along the medial axis path 
is redefined to exactly trace the boundary oi A, B, and C, then this ordering 
must still hold. Therefore, the intersection points between the edge of disc A 
and disc C must be contained in disc B. 

Second, the radius of disc B cannot be larger than the distance between 
point b and the farthest point on the edge of disc A. Otherwise all of disc A 
will be inside disc B, making disc A not a maximal disc. 

Third, the points on the edge of disc A not contained in disc B must 
include points not part of the edge of the asymmetric lens, or else disc A will 
be completely interior to disc B and disc C and not a maximal disc. 

The remainder of the argument reduces to the following. Given a point b 
in Quadrant 4, what points on the edge of disc A are farther from point b than 
Ifar^ If a line is drawn connecting point b to point a then the line intersects the 
edge of disc A in two locations, one in Quadrant 1 or 2, the other in Quadrant 
3 or 4. The intersected edge point in Quadrant 1 or 2 is the farthest point 
on the edge of disc A to b. Tracing around the edge of disc A to the other 
intersection point, each point encountered is closer to point b than the last. 
Consider the region above a line intersecting I far and point a in Quadrant 4. 
If point b was in this region, then only points on the intersection lens edge of 
disc A are farther from point b than I far This violates the third rule above, 
so this region cannot contain b. If point b is restricted to the remaining region 
of Quadrant 4, then the distance to I far is alway greater than any other point 
on the lens. Therefore, it is the case that the entire intersection overlap region 
of disc A and disc C is contained in disc B. 

The following theorem immediately derives from Theorem [3j 

Theorem 4 Let be a set of maximal discs on M(G). Let M{G) be di- 
vided into two loci of connected points Pi and P2 such that the only points Pi 
and P2 have in common are a finite set of points occupied by maximal discs 
Rn, boundary C Rn ■ L^t Rn,i be the maximal discs whose centers are on Pi 
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but do not include Rn .boundary Likewise, Let Rn.2 be the maximal discs whose 
centers are on P2 but do not include Rn, boundary The area covered only by 
the set of discs in Rna is the same for all possible Rn.2, O'^d the area covered 
by only by the set of discs in Rn,i is the same for all possible Rn.2- 

This theorem is illustrated in Figure [6j 




Fig. 5 A construction to show that disc B must contain all of the overlap between disc A 
and disc C. 




Fig. 6 If disc A is kept fixed in position, then it divides M(G) from Fig.|4jinto two parts, 
indicated as solid and dashed lines. There is no intersection between discs on the two parts 
of M{G) that is not covered completely by disc A per Theorem|3] The two parts, therefore, 
act as independent spaces. 
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3.1 Properties of an optimal planar filling 

We make some observations about the function (/)(i?jv, G), where the discs of 
i?jv have centers at points Xi e M(G). 

Per Eqn. [2] is a function of the area of discs and the area of overlap 
between discs. As the radius function is continuous over M[G) and the area of 
overlap between shapes is continuous with respect to inflating or translating 
the shape, </) is a continuous function. 

The change in (j) due to moving a single disc center is equal to the change 
in the area uniquely covered by the disc. This uniquely covered area can be 
expressed as the area of the disc minus the overlap between the disc and its 
neighbors, 0„(a;,r), where x is the position of the center of the disc D and r 
is the radius of the disc. If the disc center is moved along a path parametrized 
by i, then 

d<t> _ dMv{Vi{D)) _ 1 / dr f dOn{x,r) dx dOn{x,r) dr\\ 
dt ^ dt ~A^\ ''''di ~ dt ]h dt ) ) ' 

(3) 

The function ^ is discontinuous when ^'^'glf , or ^'^"g^'^'' is 

discontinuous. The radius function can be first-order discontinuous at a fi- 
nite number of points, (e.g. branch points) as can x(t) (e.g. branch points or 
any point of infinite curvature). At these points (j}{t) can also be first-order 
discontinuous. 

If a point of first-order discontinuity is also a local maximum (i.e. ^ 
changes in sign at the point), then small displacements of the neighbors may 
shift the sided values of the discontinuity without affecting the sign change or 
shifting the position of the maxima. The point of first-order discontinuity cre- 
ates a center trap. The local maximum at a center trap tends to stay stationary 
unless there are large rearrangements of the points in its neighborhood. As a 
result, unlike other points on M{G), a center trap tends to be a commonly 
occupied point in locally maximal filling solutions, optimal filling solutions, 
and even a fixed feature of filling solutions when N is large. 

A point of infinite curvature in M{G) that is not a branch point, (i.e. its 
maximal disc is not in contact with the boundary three or more separate but 
contiguous sets of points) is an end-point with connectivity two. 

Another point of first-order discontinuity is the point where one disc first 
contacts another. These are the points where ^'^"(^■'') or jg discon- 

^ ox or 

tinuous. Here the overlap area with another disc changes continuously from 
zero to positive. As a sign change cannot occur at this point, it cannot be an 
isolated local maximum and does not act as a trap. 
We now introduce a new term, junction point. 

Definition 5 A junction point is a point in M{G) where, relative to some 
path in M{G) that includes the point, either the radius function or the path 
itself is first-order discontinuous, or both. Junctions can act as center traps. 
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Therefore, is a continuous, piece- wise first-order continuous function. 
Insights into the structure of 4> and M{G) permit the design of an efficient 
heuristic with a high hkelihood for finding an optimal filling solution. 



3.2 Polygons 

For a convex polygon, M{G) is composed of only line segments. For a simple 
polygon, M{G) is composed of line segments and parabolic curves (Figure [t]). 
For a convex polygon, the parents are always two straight edges. For a simple 
polygon, parents can be two straight edges, a straight edge and a reflex point, 
or two reflex points. The resultant branches of M{G) and corresponding radius 
functions are given by the following three cases. 



Fig. 7 Shown are the M(G) of two simple polygons (green online). The dots (red online) 
represent junction points. 



— Case 1 (Figure ^^a)): two straight edges A line segment with a linear (or 
constant) radius function. The path and radius function can be parame- 
terized as {x{t),y{t)) = At + B, and r{t) = ct + tq, for t > 0, for constants 
c, ro > 0. 

— Case 2 Figure ^b)): a straight edge and a reflex point A parabolic curve 
with a non-linear radius function. The path and the radius function can be 
parameterized as {x{t),y{t)) = {2rot,rot'^) and r{t) = ro(t^ -I- 1), where tq 
is the minimum of the radius function. The curvature can be parameterized 
asK(0 = (l/2ro)(l + t2)"'/'. 

— Case 3 (Figure^c)): two reflex points A line segment with a non-linear 
radius function. The path can again be parameterized as {x{t),y{t)) — 
At + B. If < = is the point on the medial axis halfway between the two 
reflex points, then r{t) = ^ d?- -I- (a^jt^ — + ("y^^ ~ ^y)^^ where a is 
the distance from the halfway point to the reflex point, A = (ax, fly) and 

B=(6„6y). 
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Parent 1 




Parent 2 



(a) 



Parent 1^f^r(t) = ro(1+t2v, 



1 ro 




ir—^ 


d(t) = ro(t(1+t2)i/2 


1 


+ sinh-''{t)) 






Parent 2 




(b) 





Parent 1 




r = (ro2+d2)i/2 
/''^ ^ Parent 2 



(c) 

Fig. 8 The three types of branches (green) in a polygon are determined by the parents of 
the branch. The parents are (a) two edges (Case 1), (b) an edge and a reflex point (Case 
2), and (c) two reflex points (Case 3). For each branch, the minimum value of the radius 
function is a and the point is represented with a closed dot. The open circle represents a 
disc center at a distance d along the branch. The radius function in terms of ro, d, or a 
parametrization is provided. 
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For simple polygons, the point where a line segment meets a parabolic 
curve (or a parabolic curve joins another parabolic curve) is not a branch 
point, but does involve a change in the geometry of M{G) [7]. In Figures |8]^a) , 
[Sj^b), andjsjjc), the three cases are illustrated with the parents labeled and the 
radius function shown as a function of the path or segment length, d, or of 
parameter t. Only branch points in simple polygons are junction points. 



3.3 Center-occupied junction points and optimal solutions 




(a) (b) (c) 

Fig. 9 In (a) the labeled diagram of the isolated medial axis structure created by three 
connected polygon edges and a fixed disc is shown. For 9i = 02 = 27r/3 and t = 0.2, two 
locally maximal solutions, (b) a symmetrical solution with a disc on the L and R branch and 
(c) an occupied junction LJ solution are shown. The second solution is the global maximum. 



The strategy of the optimization algorithms in the subsequent sections is 
to find the optimal filling solution, or the global maximum of 0, by generating 
a population of local maxima. Junction points are special points that need to 
be handled separately from branches by a maximum finding search method. 

To demonstrate the diverse landscape of these maxima, and the role of 
junction points therein, we employ the following numerical experiment. A sec- 
tion of a polygon G is constructed from a ray, a segment, and another ray 
having two internal angles, 9i and 62- The medial axis of this shape is two 
straight branches (parents are a ray and the segment) and a straight branch 
that is a ray (parents are the two rays) that meet at a junction point, as shown 
in Figure Igja). The ray branch is parameterized by defining t = to be the 
junction point and defining f = 1 to be the point on the ray where a maximal 
disc would be tangent to a disc occupying the junction point. We fix a disc 
at the t position on the ray branch for < t < 1 (e.g. the large red discs in 
Figure [ojb) andjoj^c) with the x at their centers.). This region of the polygon 
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la MM 

lb MJ 

3a MJ>LJ=RJ 

3b LJ=RJ>MJ 

4 LR > LJ=RJ > MJ 

5 LR(asym) > LJ=RJ > MJ 
6a LR > LJ=RJ > LL=RR > MJ 
6b LR > LL=RR > LJ=RJ > MJ 



Fig. 10 Case (A): For the constructed problem of Figure|9] d = d\ = 62, on the top left is 
shown the form of the global maximum as a function of 9 and t. On the right the topological 
type of the global maximum is shown. The fixed disc is drawn with a dashed circumference. 
The two added discs are solid. On the bottom left, the number and form of the maxima in 
each part of the phase diagram is indicated. A key for understanding this diagram is on the 
bottom right 



below the disc is now completely isolated from the rest of the polygon per 
Theorem m We now search for all the local maxima that can be constructed 
by adding two discs below the fixed disc. We perform this search for two cases, 
(A) angles 61 and 62 between each ray and the segment are set identically to 
6* for 6* e [0,7r], and (B) 6*1 = tt/2 and 6*2 = 6* for 6* G [0,7r]. In Figure |t];b) 
and[9]jc), two examples of local maxima are shown for a case (A) construction 
where 9i ^ 62 = 27r/3 and t — 0.2. Figure [9]j a) shows the graph structure 
defined by the medial axis of this construction, with branches labeled L, R, 
and M, connected by junction point J. All local maxima of this construct can 
be classified by a two letter code representing the branches or the junction 
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Fig. 11 Case (B): Diagrams are shown for the constructed problem of Figure |9| but with 
01 = tt/2, 62 = 9- On the left is shown the form of the global maximum as a function of 
and t. On the right the topological type of the global maximum is shown. The fixed disc is 
drawn with a dashed circumference. The two added discs are solid. 



point that the disc centers occupy. We shall refer to the two letter code as the 
form of the maximum solution. Figure 9lb) is a LR maximum (one disc on the 



L branch, the other on R) and figure^ 



'c) is a LJ maximum (one disc on the L 



branch, and another on the J junction point.). In this case the LJ maximum 
is the global maximum. A RJ maximum would be symmetrically equivalent, 
and the global maximum is degenerate. 

In the top left of Fig. [lO] a diagram showing the forms of the global max- 
imum are shown for case (A) for 6 — 0.0 to tt radians and t = to 1. In the 
region labeled (LJ=RJ), the global maximum is degenerate, as in figure [9](c) . 
We observe that an occupied junction point is part of the global maximum 
for a large region of the phase diagram (LJ=RJ and MJ). In the top right 
of Figure 10 the topological types of the different global maxima regions are 
shown, each region shaded a different color. The dotted circles in the diagram 
represent the fixed disc. The solid circles represent the other two added circles. 
Two global maxima are considered topologically equivalent if they have the 
same type of edge intersections. The MM and MJ regions, for example, have 
the same type of topological intersection between the three discs. The LJ=RJ 
region has two types (colored with two shades of pink), and the LR region has 
five (colored with five shades of orange). 

In the LR region of the diagram, the form of the global maximum is sym- 
metric (similar to figure |9]^b)), except at high t. In this region, the disc on the 
L (R) branch is larger and overlaps the fixed disc, while the disc on the R 
(L) branch is smaller and overlaps only the other non-fixed disc. This region, 
therefore, also has a degenerate global maximum form. 
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In the bottom left of Figure 10 the number and forms of local maxima 
found in each part of the diagram is shown. The number of local maxima 
ranges from one to six. On the bottom right is a key for indicating which 
forms of local maxima can be found in each region. For example, in region 6b, 
there are six local maxima. The symmetrical LR maximum has a higher filling 
than a LL maximum, which is symmetrical identical to an RR maximum. The 
LL maximum has a higher filling than an LJ maximum which is identical to an 
RJ maximum. An MJ maximum has the lowest filling of all the local maxima 
in this region. Notably, the MJ local maximum is present everywhere on the 
phase diagram, except in region la, where MM is the single only maximum. 
As 9 and t both approach zero, more local maxima emerge. 

On the left of figure 11 the form of the global maximum is shown for case 
(B) where 9i — tt/2 and 02 = = 0.0 to tt radians and t = to 1. For case (B), 
there is no symmetrical solution region (aside from a region of zero area where 
6*2 = 7r/2 The space is divided into eight forms of global maxima. The only 
global maximum solution form that is not found is the RM solution form. On 
the right of Figure 11 the topological type of the global maximum is shown. 
The solution regions MM, MJ, RR, and LL each have a single topological type, 
the solution regions LJ, RM, and RJ have two, and the solution region LJ has 
three. Interestingly, the case (B) LR region has two fewer topological forms 
than the case (A) LR region from figure 10 For case (B), a diagram of the 
number of local maxima as a function of 9 and t is not shown, as we found 
the landscape too complicated to map. 

The number of maxima and solution forms and topological types depicted 
in both Figure 10 and 11 are surprisingly complicated. While the boundaries 
between different solution regions may correspond to some analytical expres- 
sion, the expression is not known. If the number of possible forms is small, it is 
simplest to check each form to find the global maximum. We observe that the 
junction point on M{G) is occupied in the global maximum solution for a large 
fraction of the diagram. This numerical experiment justifies the treatment of 
junction points as special points in the medial axis point set. This numerical 
experiment also discourages searching for an analytical expression to finding 
optimal filling solutions for at least small finite sets of discs. 



4 Algorithms for Generating Filling Solutions 

We know of no analytical method for finding the optimal filling solution for 
a given N. Instead we introduce two algorithms that explore the objective 
function landscape of to search for the global maximum. 



4.1 Genetic Algorithm 



A genetic algorithm (GA) [10] is employed to find the optimal filling solution 
for polygons. The benefit of the GA is that it uses a minimal number of 
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mathematical assumptions about the space of the fiUing solutions, although 
the computational time required is prohibitively long for N > 20. 

4-. 1.1 Algorithm Description 

GAs start with an initial random population of solutions that are combined 
and mutated until no better solutions are found after a fixed number of itera- 
tions. 

For this implementation, lOOA^ to 400A^ population members are initial- 
ized. Each member is an ordered set of coordinates of N discs. If a disc is 
randomly generated outside the polygon, it is moved inside. 

Without loss of generality, but with dramatic improvement of the efficiency 
and accuracy, the GA assumes that solutions will consist of maximal discs and 
attempts to construct them. First, the radius of each disc is grown to touch 
the nearest edge in the polygon. Second, if it can be determined that the disc 
is in a corner of the polygon (e.g. the nearby medial axis is a straight branch 
terminating in an end point), then the disc is moved to the nearest point on 
the medial axis, constructed by generating the bisector of the corner's internal 
angle. 

These constraints are applied to the entire population and then (f> is com- 
puted for each member. The population is then sorted by (j) to produce a list 
of ranked solutions from best to worst. Members of generation g are randomly 
chosen as parents for generation g + I. The relative probability p of a given 
member being chosen is weighted by its rank r, p — l/y/r. The next generation 
of members is created from the current generation as follows. 

— Best The best members, unmodified, are included. 

— Mutation One "parent" member is randomly chosen and is randomly mu- 
tated by either moving a disc randomly, displacing a disc up to 1/2 the 
polygon's width, w, displacing a disc up to l/200th of w, or moving a disc 
to a junction point point. The mutated "child" member is included. 

— Crossover Two parent members are selected and their discs are spatially 
sorted hy x*w + y, where {x, y) is the position of a disc. A crossover point 
C G [l,iV] is randomly determined. The child member contains the discs 
with indices from 1 to C from the first parent and the discs from C -f- 1 to 
N from the second. The child member is then included. 

The fraction of the next generation created by methods above can be mod- 
ified to improve the outcome of the algorithm. 

Even with the GA's capability of exploring many local maxima to find 
the global, it can still get trapped in a local maximum. The GA is run with 
different random number seeds ten or more times for each shape and value of 
N. The best of the best solutions obtained over all these runs is selected as 
the GA's final answer. 

The mathematical assumptions the GA uses are first, per Theorem [l] that 
solutions should consist of maximal discs. Second, the GA also randomly places 
discs centers on junctions. As discussed in sections |3.1| and |3.3[ junction points 
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can act as center traps and are occupied as part of many local maxima. How- 
ever, the basin of attraction around the junction point can be small enough 
that small displacement mutations do not find it. 

4.2 Heuristic Algorithm for Filling a Polygon 

In this section we introduce a heuristic algorithm that generates a putatively 
optimal filling solution of N balls, by exploiting the properties of the M{G) 
structure to generate a collection of unique local maxima. If enough are gen- 
erated, the global maximum is among them. For the two-dimensional filling 
problem, we propose a local maxima generating strategy whereby centers are 
distributed onto M{G) and the local filling maxima for that initial guess is 
found by simple gradient methods (e.g. active set or sequential quadratic pro- 
gramming optimization schemes). We also propose a method for reducing the 
number of such distributions needed to find an optimal solution by using the 
TV — 1 filling solution to generate the N filling solution. 

The first step is to intelligently divide up M{G). The medial axis is divided 
into K pieces, maximally long branch sections with monotonically increasing 
radius functions and the junction points connecting them. To generate these 
pieces for a polygon. Case (2) and Case (3) may need to be divided into 
separate sections and joined with other branch sections. The medial axes of 
the left and right polygons depicted in Figure [7] are composed of seven and 
seventeen pieces, respectively. 

Definition 6 A way, W, is a distribution of N discs over the K pieces (branch 
sections and junctions), W — {rii}^ where N — and rii E N. If the i-th 

piece is a junction point then Ui S {0, 1}. 

Conjecture 1 There is at most one local maximum per way. 

If Conjecture [T] holds, then to find the optimal N filling solution, a max- 
imum must be generated for every way of N discs and K pieces. If J is the 
number of pieces that are junctions, then the number of maxima to be searched 
is of order O {N^-''-^) (see A-4). 

Conjecture 2 Given the optimal way of distributing A'^ — 1 discs, {n'^}^ , the 
optimal way of distributing N discs is nearby, where nearby means | 
ni — n[ I is small, and that if the discs assigned to a given piece is decreased, the 
pieces that have discs increased have a minimal distance (counted by number 
of connecting pieces) to the decreased piece. 

Conjecture 3 For a given G and M(G), there is an A''' such that for A^ > A'', 
the junction points are always occupied in the optimal filling solutions. 

Given the way of the A^ — 1 filling solution, the heuristic generates the local 
maxima of the nearby ways using a local maximum finding technique. The 
best local maximum found is presumed to be the optimal A^ filling solution 
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for the shape. This heuristic is made more efBcient by taking advantage of 
center occupied junction points and the dependence of the fiUing function on 
the nearest neighbors. We implement this heuristic for polygons, which have 
a limited set of medial axis parameterized pieces to be considered. 

4-.2.1 Detailed Description of Heuristic 

Following is a more detailed description of the Heuristic Algorithm (HA). 

Auxiliary Algorithms The following sub-algorithms are needed to deploy the 
HA. 

1. Generating and Dividing M{G) into K pieces. The medial axis of the poly- 
gon is generatecj^ Parabolic curves (Case 2) and straight curves (Case 3) 
that include a minimum in the radius function are split at the minimum. 
The split branches are then recombined to form maximally long paths 
with monotonically increasing radius functions. Branch points are sepa- 
rated from branch pieces as junction point pieces. 

2. Calculating the Area of a Union of Discs. The total area of the union of 
the discs is determined analytically by dividing the space into intersec- 
tion regions defined by boundary arcs and calculating the area of each 
region [12]. The method can be applied to the entire set of discs, or, far 
more efficiently, by dividing the calculation over the discs on each piece. 
In this latter method, first the area of the union of discs on each piece is 
calculated and summed. This sum over-counts the overlaps between unions 
of discs of different pieces. Second, the overlap between the disc at the end 
of a piece and its neighboring discs on other pieces is subtracted from the 
total, once for each time it was over-counted. This latter method is more 
complicated, but also more computationally efficient because the areas of 
smaller sets of discs are calculated. 

3. Partitioning a Graph by Occupied Junctions. By taking advantage of re- 
gions of the M{G) graph isolated by occupied junctions, per Theorem |4j 
filling solutions can be divided into solutions of independent sub-spaces 
of M{G). Using the topology of the M{G) graph and a set of maximal 
discs Rni this algorithm step divides the graph into parts, or sets of pieces 
isolated from each other by occupied junctions. 

4. Finding the Local Maxima. A solution set of discs can be uniquely rep- 
resented by the way W — {ni}^ and a set of parameters {tij} where 
i e [l,ni] and j G [l,K] and tij G [0, l]^ Given an initial guess way 
and a parameter set, an optimization method (e.g. active set or sequential 
quadratic programming optimization schemes) is applied by using (jj as the 
objective function. If the initial guess includes a trial disc insertion into a 

^ using the matlab software package MatlabMedialAxis- Version 2.0 provided by Suresh 
Krishnan 1111 

^ if a piece terminates in a junction point at t=0 or 1 or both, then tij £ {0,1], tij £ [0,1), 
or tij £ (0, 1), respectively. 
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piece of a given part of M{G), then all Uj parameters of the part are free 
parameters. All Uj parameters outside the part are fixed. 

Heuristic Algorithm Given the — 1 solution; 

1. For each part of the M{G) graph isolated by occupied junctions, a new 
disc is inserted into each piece and the best solution for the part is found. 

2. For each occupied junction point of the — 1 solution, the disc center 
is removed from the junction point and initial guesses are generated by 
inserting two discs into nearby pieces and finding the local maximum. The 
more combinations of nearby pieces are tried, the larger a neighborhood is 
considered. 

3. The A^ filling solution is constructed from the best trial solution found. If 
a piece k has a parameter value t = or 1, indicating that the junction at 
the end of the piece has been occupied, then the disc is moved from piece k 
to the junction point piece. If a solution was generated for a part of M{G) 
and not included in the best solution, and the part is found in both the A^ 
solution and the A^ — 1 solution, then the solution is cached. 

4-2.2 Algorithmic Efficiency 



N=0 
N=1 

N=2 



(a) 

loooool 



10000 01000 00100 00010 00001 



20000 11000 01100 10100 10010 00011 
01001 00200 00002 00110 10001 00101 01010 




20000^ 11000 ortoo loToe-j^o 00011 

01(}01 00200 00002 00110 10001 OOIOIOTOIO 




N=2 



11000 O1Y00" 10lTre-~lQ010 00011 
100 00002 00110 10001 0010?~0T010 




20000 11000 01100 T0100 10010 00011 
01001 00200 00002 00110 10001 00101 01010 



Fig. 12 At the top of the figure is a medial axis with five pieces, three branches and two 
junctions, (a) The full table of ways is shown for N=0, 1, and 2. In (b) the search space 
is reduced using the greedy assumption that the next best solution is related to the last 
best solution, (c) We also add searches that deoccupy junction points and inserts discs onto 
nearby branches, (d) If the best 1-way was not searched, two of the four remaining 1-ways 
would have searched the best 2-way on the next iteration. 
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g)1413 31513 314lg) 



41513 31613 31514 



Fig. 13 Assume that the junction points 2 and 4 stay occupied. To generate the A'^=13 so- 
lution, three ways are searched for a local maximum. To generate N=1A, only one additional 
way, 31613, needs to be searched. The ways 41513 and 31514, can be created by combining 
the search of branch 1 and 5 with the solution of branch 3 for 7V=13. The occupied junction 
points isolate the solutions on each branch from solutions on the rest of the medial axis. 



Using Greediness to Reduce the Search Space The heuristic exploits Conjec- 
ture [2] to reduce the number of ways to search for local maxima (i.e . number 



of initial guesses) from O {N^-'^'^) to 0{N{K + J)). Figure 12 depicts a 



hypothetical medial axis with two junction points and three branch pieces. 
Figure 12 'a) shows a table of all possible ways for A^=0, 1, and 2. Rather than 
search each way of iV = 2, a reduced set is searched. That set is generated 
as follows, given the best — 1 way, one disc is added to each piece (that is 
not an occupied junction) as shown in Figure 12 b). Then, for each occupied 
junction point in the N —1 best way, the junction point is deoccupicd and disc 
is inserted into two of the branches nearby the junction, as shown in Figure 



12 |c). The maximum number of ways that will be searched, given the best 



found — 1 way, is K + AJ , where A is a constant dependent on how large 
a neighborhood of a junction point one chooses. While this heuristic is not 
guaranteed to find optimal solutions, it finds a putatively optimal N filling 
solution with only 0{N{K + J)) number of searches. 

For example, for a triangle with K = A, finding the best arrangement of 
TV = 10 discs means searching 121 ways, and the best arrangement of TV = 100 
discs requires searching 10,201 ways. By using a heuristic that exploits Con- 
jecture [2j to finding the best arrangement of iV = 10 discs requires searching 
only 70 ways, and for = 100 discs, only 700 ways. 



Applying optimization techniques to N' discs where N' < N The computa- 
tional effort to calculate the analytically exact area of the union of a set of 
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N discs is super-linear in TV, as also can be optimization methods of N pa- 
rameters. The exact order of the computational effort is dependent on the 
arrangements of the discs and the details and convergence rate of the op- 
timization algorithm. The greedy heuristic of section |4.2.2| not only requires 
searching fewer ways, but also mostly searches ways of N' discs where N' < N. 
This significantly improves the computational efficiency of finding a solution. 



Efficiently Sub-Dividing the Search Space The heuristic also improves effi- 
ciency by exploiting the properties of the solution space per Theorems |3] and 
|4]and the behavior of center traps as discussed in section [XT] 

When a junction point is occupied by a disc center in the N — 1 solution, 
the center is usually trapped and the phase space of centers can be divided into 
independent sub-spaces. If it is known (or guessed) that the best solution for N 
also includes a center at the junction, then the sub-parts of M{G) connected 
only by the junction point can be searched independently. Per Theorem [4j 
rearrangements of centers in one sub-part cannot affect the best arrangement 
of centers in another if they are connected only by a center-occupied junction. 

This means that searches can be performed on a subset of the N discs 
(efficient per section 4.2.2) and that solutions of independent sub-spaces of 
M{G) can be cached. Figure 13 shows how, when junction points are presumed 
to remain occupied, only one additional search of a way is needed to generate 
the next putatively optimal N filling solution. Per Conjecture|3) at sufficiently 
large TV, junction points are occupied. This implies that for large N , generating 
the optimal N filling solution from the optimal iV — 1 filling solution requires a 
search of only one additional way, reducing the complexity of the HA to 0{N) 
searches. 



4-. 2. 3 Self- Correcting 

When implementing the HA, a practical choice is made as to how large a 
neighborhood of ways N discs that are nearby the optimal — 1 way will 
be searched. There is a computational trade-off between searching only ways 
such that I Ui — n[ \= 1, in which case the optimal way may be missed, 
or such that | ni — n[ |< oo, in which case the optimal way can not be 
missed but the search space has not been reduced. 

One weakness of a method that uses the A^ — 1 solution to find the A^ 
solution is that, if the optimal A^' solution is not found, all solutions for N > N' 
may not be optimal as well. However, we observe that in most cases where the 
HA does not find the optimal A^' filling solution, by some N > N' , the HA is 
generating the optimal solution again. That is, even if the wrong solution is 
found, we observe that the solution finding method tends to self-correct at a 



higher A^. In Figure 12 ^d), for example, if the optimal A^=l way was omitted 
from the search, the optimal A^=2 way would still be searched by two of four 
alternate A^=l ways. 
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Presuming Conjecture [3] is true, we can show that for at least one common 
case, an HA searching only a small neighborhood of ways still always self- 
corrects. Consider, an M{G) with K pieces and J junctions, where all pieces 
that are sections of branches are connected by pieces that are junctions (e.g. 
G is a convex polygon) . Assume that we restrict our search to ways such that 
all the J junctions are occupied for N > J. The N — J remaining discs will be 
partitioned over the remaining K — J pieces. For this case, per Theorem|4] the 
increase in the filling measure due to adding an arbitrary number of discs to 
a piece can be solved independently. Also, adding discs sequentially and opti- 
mally to a given piece strictly increases but the change in cj) monotonically 
decreases. It follows that the best N way found with the given constraints can 
never include removing a disc from a piece. Thus the heuristic will search for 
the best N way by comparing the local filling maxima generated by adding 
one disc to each of the K — J pieces. It follows that the heuristic will always 
find the best N filling solution for which the J junctions are occupied. If Con- 
jecture [Sj is true, since for some N > N' all the junctions will be occupied 
in optimal solutions, it follows that despite having not generated the optimal 
way for all N < N' , the heuristic generates the optimal way for all N > N' . 

We propose a stronger conjecture than Conjecture [3] 

Conjecture 4 For a given G and M{G), there is an N' , such that for N > N' , 
the junction points are always occupied in all filling solutions that are local 
maxima. 

If this conjecture is true, then it would also follow for M{G) with K pieces 
and J junctions, where all pieces that are sections of branches are connected 
by pieces that are junctions, the heuristic will always self-correct and generate 
the optimal way for sufficiently large N. 

4-2.4 When the Heuristic Algorithm fails 

Even if the Conjectures [T] [2] and [3] above hold, in practice this HA may still 
fail to find the optimal solution for the following reasons. 

(1) Assuming a way has no local maximum While searching for the local 
maximum associated with a way, it is common to generate the local maximum 
of a nearby way instead (e.g a junction point becomes occupied). This leads 
to the conclusion that the way has no local maximum. However, the search 
may simply have been initiated outside the basin of attraction of the local 
maximum of the way. 

(2) Searching in too small a neighborhood As discussed above, some optimal 
N solutions require looking in a larger neighborhood of the TV — 1 solution. 
Tradeoffs that balance confidence in finding the optimal solution against the 
computational cost of searching larger neighborhoods may result in optimal 
solutions being missed. 
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HA and GA Best Way: Best Way: 
Way Match HA GA 


Best <j) : 
HA 


Convex 
Concave 


98.1% 1.9% 0% 
92.97% 3.4% 3.63% 


100% 
96.37% 



Table 1 Table 1. A comparison of the filling solutions generated by the HA and GA for 
five convex polygons and 21 concave polygons for Af=l-21. 



(3) Solutions are only as good as the optimization method applied Lastly, local 
maxima finding techniques can have trouble converging. This is not a failure 
of the Heuristic Algorithm, per se, but occasionally affects the HA solution. 
Switching which nonlinear constrained minimization optimization technique is 
being applied generally solves the problem. 



4.3 Heuristic vs. Genetic Algorithm Filling Solutions 

To assess the capability of the heuristic algorithm vs. the genetic algorithm, 
solutions were generated for A^=:l to 21 for a selection of five convex polygons 
and 21 concave polygons. The putative best solutions produced by this HA 
match well the solutions generated by the GA. Specifically, the HA almost 
always produces solutions of the same way as the GA. The gradient optimiza- 
tion technique employed by the HA is usually better at converging to a final 
set of disc positions for a given way than the GA. On rare occasions the HA 
and GA find different ways. When the HA way is better, the GA has usually 
become trapped in the wrong local maximum. When the GA way is a better 
solution, wc find that the way was outside the neighborhood that was searched 



by the HA. Examples of filling solutions are shown in Figure 14 



5 Optimally filling a polygon as N — >■ cx) 

It is instructive to examine how the optimal filling of a shape converges to 
the total volume of the shape as — >■ oo. As iV — >■ oo, the centers will be 
distributed densely in M{G) such that the change in the density of centers 
measured over small intervals of the 1-manifolds M{G) can be can be con- 
sidered a smooth continuous function in the continuum limit. As discussed in 
reference [T], the continuum limit solution can be solved exactly for simple 
polygons by analyzing the three types of branches found in simple polygons 
(Case (1), (2) and (3) of Figure |8];a),[8](b), and|8];c)). For the Case (3) type, no 
disc centered on such a curve fills any more area than what is filled by placing 
two discs at the ends of the curve. Thus in optimal solutions. Case (3) type 
curves are empty except for their ends. 

Let p{t) represent the density of centers along a parameterized path of 
Af (G), r{t) be the radius function, and n{t) be the local curvature of the path. 
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N =21 N =21 N =21 N =21 N =21 

Fig. 14 Examples of the optimal filling solutions of three convex and two concave polygons 
for N=l-21. The top row shows the medial axis of each polygon. 



where {x(t),y{t)) is the parameterization t g [ia,ib]- Given an expression for 
the unfilled area Ai along the path i of M{G) of the form, 



(4) 
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where C; is a function to be determined, we would like to determine the func- 
tion p that minimizes this area constrained by 

7V = ^ pdt. (5) 

Note that if we sum the unfilled areas Ai over all of M{G), then </> = 1 — 
^(Ai/Ac), where Aq is the area of G. This variational problem can be solved 
by forming the Lagrangian 

L[p{t)-\]= j'^" [c,{ny,r)^ + \(^dt (6) 

and taking the pointwise derivative with respect to p(t), 

This relationship is satisfied by functions p that satisfy 

- 2C,iK, /, r) + p-^C\{k, /, r) + p^X = (8) 
dp 



Solutions of the form p = ^^lil^^^L^^ satisfy this equation. 
For Case (1), where {x{t),y{t)) —At + B, r ~ rot, and ta > 0, 

C= (l-r'2)^^V(12r) (9) 

as shown in Appendix A-1. It follows that p =oc r^^/^. 

For Case (2), where {x{t),y{t)) = (2rot, roi^), r = ro(t^ + 1), ro is the 
minimum of the radius function, and K{t) = (2ro) ^ (l + t^) , 



.5/6 

or p oc r~^' ^. 



as shown in Appendix A-2. It follows that p = po ^j^^^ 

For both Case (1) and Case (2), the distribution of centers follows a power 
law with respect to the local radius function. Centers on M{G) will be dis- 
tributed more densely where the radius function is smaller. Given p = po?" 
for a = 1/3 or 5/6, pa can be determined from Equation [sj 

Po = N (^Jl' r-'^dty' (11) 
Po - N/Ro. (12) 

Ro is then a constant determined by the radius function of the branch section 
of M(G). 
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For Case (1), the distribution of centers on the medial axis path is also 
scale-free. The distribution of centers also follows a power law with respect to 
the distance from the vertex (where t — 0) the polygon. 

Equation |4] becomes, 

j\lC{K,r',r)dt^^C. (13) 

Thus, in the continuum limit the optimal filling solution converges to the 
area of the shape with an asymptotic error proportional to N~'^ for ideally 
distributed centers. We presume that all shapes that can be approximated by 
simple polygons with an increasing number of sides also converge with an iV~^ 
error term. 

If we divide M{G) into k branch sections we can predict what fraction of 
the discs (Ni/N) will be distributed over each branch i a.s N ^ oo. 

k 

A = ^^.(iVO (14) 

1 

k 

N^Y.^^ (15) 

1 

Since we have distributed our discs optimally, we can treat A^{N) as a con- 
tinuous function and thus 



Arbitrarily setting j = fc, 

dAi dAk C, Ck , , 

The fraction of discs on a given branch i is, 

^ ^ (19) 

N (Ci)l/3 + (C2)l/3 + ... + (Cfc)l/3- ^ > 

For a triangle, which is always composed of three Case (1) branches, the 
fraction of the discs on a given path can be solved analytically to be 

f. = (20) 

cot (6'i ) + 001(6*2)+ cot (03 ) 

where 9i is an internal angle of the triangle, each of which is associated with 
a branch. From equation |20[ it is clear that the optimal solution preferentially 
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populates medial axis branches associated with smaller internal angles. This 
can be observed in the optimal filling of a triangle in Figure |15[ 

An illuminating case to consider is a branch of M{G) where both the 
curvature of the path and the radius function are constant. In Appendix A- 
3, we find C = (k^t + or a constant over the branch. It immediately 
follows that p = N/T, where T is the length of the branch. As expected, 
centers are distributed evenly over the branch. For this case, C = T^C. Thus, 
via equation 19 M{G) that can be divided into branches of constant curvature 
and radius functions also have known distributions as — >■ c», and branches 
with higher curvatures will be more densely populated. 



6 Conclusion 

In this paper we investigated the new problem introduced in reference reference 
[T] of optimally filling shapes with balls of varying radius. Filling combines 
two classic mathematical problems, the packing of shapes p^fT4l[T5l[T6l[T7lfT8] 
and the covering of space _18, 19,20,21,22 . Like in packing problems, the balls 
cannot overlap the boundary of the shape, but, like in covering, the balls may 
overlap each other without penalty. This combination of constraints generates 
an interesting new problem. 

In our research, filling solutions arise from the problem of modeling anisotropic 
nanoparticles as rigid bodies composed of a sum of isotropic volume-excluding 
potentials 23,24,25,26 . Filling solutions have applications to many other ar- 
eas of optimization, including the problem of irradiating a tumor with the 
fewest number of beam shots, while controlling the beam diameter, but with- 
out damaging surrounding tissue |27j : using time-delayed sources to create 
shaped wavefronts; combining precision-placed explosives with tunable blast 
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radii; positioning proximity sensors with defined radii; cell phone and wireless 
network coverage; or any problem of ablation or deposition where one has a 
sharp impenetrable boundary and a radially tunable tool. 

We find the filling problem to be surprisingly rich. This paper describes 
the basic structure of the filling problem in arbitrary dimensions. For polygons 
we have provided a deeper description of the solution space and detailed two 
methods for finding numeric approximations of the optimal filling solutions. 
We also have shown how optimal solutions in polygons converge to simple 
analytical expressions as the number of discs approaches infinity. 

The solution space structure of a simple polygon has features that we 
expect to find in more generalized and higher dimensional shapes, namely 
first-order continuous manifolds that join at lower dimension manifolds where 
centers are trapped. We predict that higher dimensional polytopes will also 
have manifolds with scale-free solutions. 

Even in two-dimensions, many open questions remain. For example, how 
can optimal solutions be found for a generalized shape G that is not a simple 
polygon? Also, for a given shape G, is there a N' , such that for N > N' 
junction points are always occupied by centers in optimal solutions? Many of 
the potential applications of the filling problem demand solutions for three- 
dimensional shapes. It is desirable to develop practical methods for finding 
optimal solutions in higher dimensions. 
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8 Glossciry 

maximal ball A ball contained completely in a shape G that is not a proper 

subset of any other ball also contained in G. Also, a ball tangent to the surface 
of G at at least two points, that is completely contained in G. In a 2D planar 
shape, a maximal ball is a maximal disc. 

medial axis M{G), the locus of the centers of all maximal balls of G. 

radius function The radii of the maximal balls of a shape G. 

normal point A point on M{G) that is the center of a maximal disc in contact 
with the boundary S at exactly two separate but contiguous sets of points. 

end point A point on M{G) that is the center of a maximal disc in contact 
with the boundary S at exactly one contiguous sets of points. 

branch point A point on M{G) that is the center of a maximal disc in contact 
with the boundary S at three or more separate but contiguous sets of points. 

branch A set of contiguous normal points on a medial axis. 

parent of a branch The two contiguous parts of S from which the normal points 
of the branch are derived. For a simple polygon, parents can be a polygon edge 
or a reflex point. 

neighbor If a maximal disc has a disc center that can be reached by a path 
along M{G) starting at the center of maximal disc A without traversing a 
third disc center, then it is the neighbor of maximal disc A. 

center trap A point on M{G) where a first-order discontinuity coupled with 
a local maximum in (all centers fixed) creates a local maximum that is 
stationary with respect to small changes in the position of the neighboring 
discs. 

junction point A point on M{G) that can act as a center trap. For a polygons, 
branch points are junctions. Whether a generalized planar graph can have 
junction points that are not branch points we leave as an open question. 

piece A junction point or a section of a branch. 

way A distribution of N discs over the K pieces that compose M{G). 

part A connected set of pieces only connected to pieces not of the set by 
disc-center occupied junctions. 
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9 Appendix A-1: Distribution function along a medial axis branch 
with no curvature and a linear radius function 




We calculate the area between two discs along a medial axis branch gen- 
erated by two polygon edge parents as the two discs approach each other. 

Figure AjT] shows two overlapping maximal discs, of the same radius, sep- 
arated by a distance d, and one of the two lines tangent to both discs. The 
green region is the uncovered area in between the discs and the tangent line. 
As d — >■ 0, what is the area, A of the green region? 



A ~ rectangle — 2quarter circles ~\- ^lens (A.l) 



Using an acosine Taylor expansion and the square root 

We now approximate the uncovered area between two discs and the polygon 
edges for a radius function is that is not a constant, by bounding the answer 



Optimal Fillings 



33 




between an upper and lower bound. Make one of the discs of Figure A{T] larger 
by Ar — dr' , as per the Figure A[2]^a). The distance between the two centers 
is still defined as d. Put a cotangent disc of radius R at the point of tangency 
of the both discs per Figure Aj2|^b) and Figure [2j^c). If i? = r or if i? = r + dr' , 
the centers are now d^/1 — jr'^ apart. The uncovered area of Figure A. a) is 

1 d^'iVT^f 



and 



1 d^jVl-r'^y 



bound between Figure A. Sfb) and Figure A. 2rc) or between 



24 



r+dr' 



. As d — > 0, the area uncovered is 



1 d^{y/l-r'^f 



^uncovered 



24 



(A.6) 



The area is then doubled to account for the identical uncovered piece on the 
other side due to the other tangent line (i.e. polygon edge). 
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We now observe that d — ^ where p is the density of disc centers along the 
branch. To determine the total uncovered area along a branch of length T, we 
would sum all the uncovered areas that are at density p along the branch, or 

^f^{^-rfy^>f^{l^n'''^< (A,7) 







12rp3 Jo Ylrp- 



10 Appendix A-2: Distributions along the parabolic medial axis 
branch of a polygon 




Fig. A. 3 The area shaded in green is the uncovered area between two discs of the different 
radius and the polygon edge along a parabolic path. 



For concave polygons, the medial axis branch associated with discs tangent 
to the reflex point and an edge of the polygon is a parabolic curve. The reflex 
point forms the focus of the parabola and the polygon edge the directrix. If the 
ends of such a parabolic curve are occupied, all the uncovered area is between 
the directrix and the set of overlapping discs distributed along the branch. The 
parabola has both a changing radius function and changing curvature along 
the branch. 

First we note that for two discs that have centers of distance d apart, the 
uncovered area between the discs and the directerix is the same as Equation 

El 

Area^^(l-r-)^/^ (A.8) 

The local density of centers is p = where As is the arc-length between 

the two centers. However, as d — >■ 0, d ~ As so, 

Areac.^(l-r-)^/^ (A.9) 
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Using a parameterization of the parabola, where y = at^, x — 2at, r — 
at^ + a, the arc length s{t) = a (t^/l + + sinh~^t) , and the curvature = 

^(l + ^^)-^/^then, 

dr dr dt 2at t . . ^ 

-r^-r-r^ ; = A.IO 

ds dtds 2a^/^T^2 ^TTF ^ ' 

Note that r' < 1, which is a general property of r' of a medial axis. 
So, 

il-r''f/^=(-^\"\2a.. (A.11) 



Now Equation A. 9 is equal approximately to 

1 / 2roK 
2V V~ 



Area^^(^) (A.12) 



where tq is the smallest radius of the parabola, or ro = a. Or, 

1 / 1 



Area-- ^ -r ] . (A.13) 

If the parabola is defined from ta to ti,, then the total uncovered area is 



11 Appendix A-3: Constant curvature, constant radius function 

Consider the case of two identical discs separated by a branch of length d and 
curvature w k. What is the uncovered area between them? See Figure A|4] 
The uncovered area is the area swept between the two arcs tangent to the 
discs minus the area covered by the discs within the area swept, or twice the 
green area of Figure A|4]d. Let be the angle ZABC, 9c be angle ZDCE. Let 
del be the chord length between A and C. Then 9 ~ dn, d^i = ^sin (^), and 
dc = 2cos^^ (ct^^^ (t^))- ^^'^ s.T:<ia. of the green region of Figure A|4j3 is equal 
to the half disc minus the grey area, or ^{9c — sin{9c)). 





-Trr^ +r2(6'c- sin(6'cl)^. 15) 
2rd-Tir'^ + r'^{ec-s\n{9c)) (A.16) 



Expanding the third term 
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2r- (^cos-i (^sin (f )) - ^1 - ( J,sin (f ))^ ( ^.sin (f )) j (A.20) 

Substituting a Taylor series expansion for the inverse cosine and square 
root 

(-/2 - (;^Bin (f )) - 1 ( J,sin (f ))^) - (A.21) 

2-^(l-H;^-n(t))')(;^sin(t)) (A.22) 

= 2r^ (-/2 - (;|sin (f )) + § (^sin (f ))^) (A.23) 



Thus, 



Auncovered = 2rd ~ ^siu ^ j + fsiu f^" ) ) (A.24) 



Substituting a Taylor series expansion for the sine terms. 

^ , 4r //rfK\ 1 1 4 

12 12 r ^ 

So the uncovered area along the whole length of the branch is {d = -) 

where N is the total number of discs distributed over the branch section. 
Using equations (10) and (14) of the main paper, if M{G) can be broken 
into branch sections with approximately constant r and curvature, then the 
fractional distribution of the N discs over each section can be determined such 
that the fraction of discs distributed over a section k of length is 

1 \ 1/3 



fk cx Tfc [^nirk + — j (A.28) 

We observe that, in general, the density of discs is higher in regions of high 
curvature. 



12 Appendix A-4: How many ways? 

Per conjecture [T] to find the optimal N filling solution, a maximum must 
be generated for every way of N discs and K pieces. How many maxima 
searches is this? Assume that the K pieces have J junctions, J < K. For 
m € N, < m < J, there are ways to occupy the junctions, leaving 
N ~ m remaining discs to allocate over the K — J remaining pieces. A weak 
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composition is a way of partitioning an integer into a sequence of non-negative 
integers, where order matters. The number of weak compositions of TV — m 
discs over K-J pieces is 

Thus the number of ways to be searched is equal to X]m=o'''^^ (m) 



(A.29) 



min(J,N) / -.\ t - 

fJ\/N-m + K-J-l 
h Wl K-J-l 

' Jl{N^m + K-J-iy. 

£^^ml{J~my.{K-J-l)\iN-m)l ^ ' 

= 0{N^'-'-^) (A.31) 
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